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ABSTRACT 

We propose a computer-efficient and accurate method of estimation of spa- 
tially correlated errors in astrometric positions, parallaxes and proper motions 
obtained by space and ground-based astrometry missions. In our method, the 
simulated observational equations are set up and solved for the coefficients of 
scalar and vector spherical harmonics representing the output errors, rather than 
for individual objects in the output catalog. Both accidental and systematic 
correlated errors of astrometric parameters can be accurately estimated. The 
method is demonstrated on the example of the JMAPS mission, but can be used 
for other projects of space astrometry, such as SIM or JASMINE. 

Subject headings: astrometry — parallaxes — proper motions — methods: nu- 
merical 



1. Introduction 



Projects of global space as trometry, su ch as Hipparc os (jESAl 119971 ). Gaia (jPerryman 



2003|), SIM flUnwin et al.ll2008l ) or JMAPS flGaumell201ll ). result in very large systems of 
observational equations, which contain the unknown parameters of interest, as well as an even 
larger number of instrument calibration and satellite attitude parameters. The structure of 
this system depends on the architecture of the mission, adopted schedule of observations 
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and properties of the instrument. In its most general matrix form, a linearized observational 
equation of space astrometry is 

Aa + Cc + Pp = y, (1) 

where the three main types of unknowns are separated, namely, the astrometric parameters 
a, the calibration parameters c and the attitude parameters p. The measurements collected 
during the on-orbit operation are involved in the right-hand side y of the observational 
equations. These measurements include accidental and systematic errors. The accidental 
error arises from a completely stochastic, unpredictable process and occurs spontaneously 
in each observation, independent of the previous or subsequent observations, for example, 
photon shot noise in the incoming light signal. The systematic error has a certain determin- 
istic cause, such as a particular state of the instrument, and may in principle be avoided 
or mitigated by setting adequate requirements to the mission design and the instrument. 
The problem we are going to consider is how to calculate the propagation of accidental and 
systematic errors in the astrometric solution in the general case. In particular, we are con- 
cerned about the occurrence of correlated and zonal errors in the resulting proper motions 
and parallaxes. 

We are considering in this paper the position-correlated part of the total absolute error 
of an astrometric catalog or reference frame. It can be visualized as a smooth pattern in the 
distribution of the absolute "observed minus true" error on the celestial sphere. The actual 
distribution of absolute error is sampled on a discrete set of objects (e.g., stars and quasars), 
so that the uncorrelated part will look like discrete noise, but if a smoothing procedure is 
applied to this sampled function, the underlying large-scale pattern will emerge. This smooth 
pattern will be stochastic, i.e., unpredictable, if it is caused by the accidental errors of the 
measurements, or deterministic, and therefore predictable, if it is caused by a systematic 
measurement error. In reality, any catalog produced by a space mission will include a 
mixture of accidental and systematic errors. The goal of a well-designed astrometric mission 
is to keep the systematic error smaller that the accidental error. However, an external 
catalog of superior accuracy is usually required to verify that this goal is achieved. For 
example, a comparison of the Hipparcos catalog with the FK5 system of 1535 reference stars 
revea led zonal errors in the latt er that were larger than the estimated accidental error per 



star (IMignard fc Frceschlell2000l ). The type of large-scale errors considered in this paper is 
different from the accumulated small-scale errors, which originate in the shared attitude error 
for stars observed w ithin the same field of view (or the same scan for a scanning mission) 



(Ivan Leeuwenl 120071 ). The details of such accumulation of small-scale error in the in-scan 
direction for scanning missions is highly dependent on specific mission implementation and 
must be modeled separately on a case-by-case basis. 
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The smoothness of spatially correlated large-scale errors invokes the use of spherical or- 
thogonal functions for their representation and analysis. A parallax error is a scalar function 
of spherical coordinates, and is best described by the scalar spherical harmonics. Proper 
motion and position errors are tangential vector fields on the unit sphere, and should be 
properly represented by vector spherical harmonics. The spherical harmonics are basis func- 
tions, therefore, any smooth function of coordinates can be accurately decomposed into a 
linear combination of a sufficient number of spherical harmonics. Some of the spherical 
harmonic terms correspond to specific physical parameters. For example, the zero-order 
scalar spherical harmonic of parallax is a constant function (unity) on the sky. This term 
is the well-known parallax zero-point. The zero-point is the average of all absolute parallax 
errors of a given catalog, which is different from zero. This common type of error is of spe- 
cial significance in astrometry, because it directly affects the cosmic distance scale based on 
trigonometric parallaxes. A common offset of parallaxes is relatively more important for dis- 
tant objects, hence, the empirical calibration of such mostly distant calibrators as Cepheids 
or RR Lyr-type stars is especially vulnerable to a zero-point error. All higher-order harmon- 
ics do average to zero on the sphere and have more localized effect. As another example of 
coupling between physical parameters and specific large-scale errors, the Oort's parameters 
describing the differential rotation of the Galaxy are represented by certain vector spherical 
harmonics of the proper motion field (IMakarov fc Murphyl 120071 ) . 



One might naively assume that random uncorrelated noise in the observational data 
results in a uniform distribution of error among the different types of zonal error, i.e., in a 
flat spherical harmonic spectrum. In fact, this is never the CcLS6, clS the mode of operation, 
the geometrical order of observation, and the finite field of view of the telescope, all lead to 
a strongly non-uniform propagation of observational noise in different orders of correlated 
error. The technique presented in this paper goes back to the idea of using orthogonal 
functions for the analysis of on e-dimensio nal "abscissae" errors of the Hipparcos mission 
(for definition of abscissae, see lESAl Il997l . Vol. 2). The concept of Hipparcos was built 
around a stab le and self-calibrat ing basic angle separating the two viewing directions of 
the telescope. iHoyer et al.l ( )198ll ) suggested that the problem of propagating perturbations 
of the basic angle beco mes tract a ble if these perturbations are represented by a Fourier 
series. Using this idea, iMakarovl ( 119921 ) showed that Fourier time-harmonics of the basic 
angle propagate with uneven magnification coefficients into the corresponding harmonics 
of star abscissae, with the 60 deg-harmonic being the dominating one in the output. The 
reason for this peculiar propagation of white noise is the value of the basic angle, which is 
close to 60 deg. The technique was further developed and applied to the R0mer project^] 



lr The R0mer concept was the precursor of the Gaia mission. 
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by iMakarov et al.l (119951 ) . The propagation of periodic perturbations of star abscissae into 
specific spherical harmonics of absolute error is a more complicated issue, which depends 
on the scanning law and the geometry of r eference great c i rcles. This generalization was 
implemented for a Hipparcos-like design by IMakarov et al.l ( 119981 ). where the 2D errors of 
parallax were represented by spherical harmonics. Using an orthogonal basis of functions on 
the unit sphere puts the concept of absolute astrometry on a rigorous mathematical footing. 
The proposed JM APS and SIM missions, u nlike Hipparcos and GAIA, employ only one 
viewing direction ( jZacharias fc Dorlandl 120061 ). The emergence of correlated errors in these 
projects is defined by the density of overlapping observational frames and the accuracy of 
quasar constraints. 



2. Observational equations and quasar constraints 

Linearized observational equations of global astrometry relate astrometric parameters of 
interest, satellite attitude parameters and instrument calibration parameters to astrometric 
observable parameters. Although the actual relations between these parameters are nonlin- 
ear, the required linearization can be achieved by taking the perturbation form and limiting 
the Taylor expansions on the left-hand side of the equations to first order. If the initial guess 
or prior knowledge of the fitting parameters is close to the truth, this first-order approxima- 
tion provides an accurate solution; otherwise, the linearization and global solutions have to 
be iterated. 

A common property of all proposed space astrometry projects is that most of the nui- 
sance parameters (e.g., calibration and attitude) entering the observational equations should 
be determined from the same observations along with the star parameters. Essentially, an 
instrument for space astrometry is self-calibrating and self-navigating. Some of the crucial 
parameters, e.g., the basic angle for Hipparcos, the baseline length for SIM, are stable by 
engineering requirements; they should be re-determined relatively infrequently during the 
mission. On the contrary, the attitude parameters are unique for each astrometric frame 
or scan, and therefore, generate the bulk of nuisance parameters. The Euler angles (or 
quaternions) of spacecraft attitude are approximately known from the navigation system, 
including a separate star tracker device, but much more accurate values of these parameters 



are determined from the main observations themselves (ILim et al.l 120101 ) . 



It was shown on the example of the SIM project that coupling between the attitude 
unknowns and the astrometric unknowns can cause a loss of condition and a non-uniform 



propagation of errors in a global solution ( IMakarov fc Milmanl 120051 ) . A strict relation exists 



between the basis vectors of parallax distribution, obtained by the singular value decom- 
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position (SVD) of the corresponding part of condition equations, and the scalar harmonics 
sampled on a discrete set of stars. The reciprocal singular values are simply the coefficients of 
different degrees and orders of error, propagating into the final parallax solution. By virtue 
of the relatively small size of the design matrix, the SIM grid solution was ideally suited for 
rigorous mathematical analysis of various aspects of error propagation. Other astrometric 
missions invoke much larger Least Squares (LS) problems, and SVD analysis becomes in- 
tractable. In this paper, we are setting out to develop a numerical method to estimate the 
propagation of large-scale correlated errors in very large LS solutions. 

In the perturbation form, the unknowns in Eq. [1] are small corrections to a priori 
parameters describing the stars and the state of the instrument, and y is the vector of small 
differences between the the predicted and the actual measurements. The grand design matrix 
can be constructed, in principle, from the individual blocks A, C and P, although it is never 
done in practice because of its huge size. Th e standard method of solving such problems 



is block adjustment (e.g. Ivon der Heidd 119771 ). using the natural sparsity and structure of 



the design matrix. Briefly, there have been two algorithms considered for large astrometric 



probl ems, the iterative block adjustment and the global direct solution (iBucciarelh et al. 



199ll ). Hipparcos, the only implemented astrometric space mission thus far, relied on the 
iterative adjustment, in which the major blocks of unknown parameters were estimated and 
updated in turns, while keeping the other types of parameters fixed, resulting in a number 
of iterations across the range of mission parameters. The convergence of the large-scale 
iterations can not be taken for granted, but should be verified by dedicated simulations. A 
similar algorithm of iterative block adjustment has been developed for the Gaia data analysis 



system (ILammers et al.ll2010l ). On the mathematical grounds, a global adjustment, which is 
a simultaneous, one-step solution for the multitude of mission parameters, should be more 
exact, faster and easier to analyze, but it poses a considerable implementation challenge for 
huge LS problems with a large number of nonzeros. 

The block structure of the grand design matrix is defined by two different types of 
dependancies of the unknowns. The astrometric unknowns a are object-dependent, i.e., they 
comprise independent sets of several unknowns for each object. Five astrometric unknowns 
per star are usually considered, namely, position components, parallax and proper motion 
components. The calibration unknowns c and the attitude p are mostly, but not exclusively, 
time- dependent. If the observational equations are sorted by time, the nonzero elements of 
the design matrix corresponding to the attitude unknowns are found in tight, relatively small 
blocks, because the attitude of an astrometric instrument is fast-varying. The calibration 
parameters can be discretized too, so that a separate set of calibration unknowns is assigned 
to a fixed interval of observations, which we will call a calibration block in this paper. The 
calibration parameters are expected to change slowly with time. Therefore, the corresponding 
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blocks of nonzero elements are longer that the attitude blocks. It is convenient to adjust 
the discretization steps in such a way that the boundaries of the blocks are aligned, so that 
an integer number of attitude blocks corresponds to each calibration block. The astrometric 
unknowns in such a design matrix are scattered across its entire length, because the same 
object is observed multiple times during the mission. 

It is sometimes practical to eliminate the attitude and calibration unknowns in the 
equations rather than solve for them directly along with the astrometric unknowns. This 
elimination is achieved by the QR factorization of each block and the subsequent Q r or- 
thogonal transformation of the remainder of the design matrix and the right hand side of 



the equations, as described in (IMakarov k, Milmanl 120051 ). Because of the nested structure 



of the blocks, the elimination procedure becomes hierarchical, the smaller attitude blocks 
eliminated first, followed by the calibration blocks. As a result, the number of unknowns is 
significantly reduced. However, this reduction comes at a cost, because the design matrix 
becomes much denser. Obviously, nonzero off-diagonal elements are generated for any pair of 
objects, which were observed within the same calibration block. The degree of densification 
depends on the average number of objects within a calibration block. To avoid intractably 
dense matrices, smaller calibration partitions are preferred. In the JMAPS global solution, 
several large-scale calibrations parameters are solved for each frame, along with the three 
attitude unknowns. In that case, the direct LS solution is obtained for about 29 million 
unknowns with 1 million grid objects, or 34 million unknowns with 2 million grid objects. 
After the QR elimination, only 5 or 10 million unknowns remain, respectively, but the design 
matrix is much denser. The number of equations to be solved is 144 or 288 million, respec- 
tively. Our idea presented in this paper is that in many cases, it is sufficient to consider the 
correlated errors of the simulated mission, rather than the individual errors of numerous grid 
objects. This strategy helps to reduce the number of astrometric unknowns to manageable 
levels, fully capturing an important characteristic of mission performance. The mathematical 
technique is described in § §2.11 and 12.21 



2.1. Spherical harmonics 

The astrometric part of observational equations can be written as 

(5s + 8fi(t - t ) + M( s o • b) s - b)) • r = y. (2) 

where 8sq and 5fi are the unknown corrections to mean position and mean proper motion at 
t Q , tangential to the celestial sphere at s , which is the assumed position unit vector of the 
objects at time to, bw is the unknown correction to parallax, and b is the position vector of 
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the space craft with respect to the barycenter of the Solar system at the time of observation t, 
assumed to be known. The vector r is a certain fiducial direction defined by the instrument, 
for example, the baseline vector of SIM, or the nominal scanning direction of Gaia. This 
vector depends on the instantaneous attitude and the calibration parameters, but here it 
is assumed to be known, because all the nuisance parameters have been separated in the 
linearized equations into independent blocks. For a two-dimensional pointing mission such 
as JMAPS, two condition equations emerge from a single observation, since there are two 
fiducial directions in the focal plane, corresponding to the rows and columns of the detector 
array. The right-hand side of Eq. [2] includes the measurement and additive accidental and 
systematic errors. 

The astrometric condition equations are linear and can be solved by direct LS with or 
without elimination of the nuisance parameters. The main technical problem arises from the 
size of the normal matrix, which require supercomputing facilities and advanced algorithms. 
For Gaia, the size is so large that a dire ct LS solution is deeme d impossible, and the adopted 



iterative solution still takes a long time ( lO'Mullane et al.ll201ll ). Solving directly for up to 34 



million unknowns has been proven feasible with a specially adapted PARDISO solver (PAR- 
allel Direct SOlver, part of the Intel Math Kernel Library), but it still takes several hours 
of computing time to complete. For the testing and verification purposes, full-scale runs of 
the global solution have to be performed multiple times, with various input data. Our idea 
presented in this paper is that in many cases, it is sufficient to consider the distribution of 
error on the sphere, rather than individual errors of numerous grid objects. Thus, we substi- 
tute the object-dependent astrometric unknowns in Eq. [2] with the expansions in spherical 
harmonics, which are functions of celestial coordinates, e.g., the ecliptic coordinates A and 
0: 



oo 



5s = 'Sj V(A, (3) position 

3=1 
oo 

= MjVj(\, (3) proper motion 
i=i 

oo 

5vjq = Qj Uj(X, 0) parallax (3) 
with Uj being the scalar spherical harmonics and V, the vector spherical harm onics. For 



a detailed description of spherical harmonics see, e.g., ( iMakarov fc Murphyl 120071 ). Here we 
only reproduce some basic formulae. The vector harmonics Vj are composed of two types 
of functions, called magnetic and electric harmonics, and E™ respectively. These vector 
harmonics are derived via partial derivatives of the scalar spherical harmonics over angular 
coordinates, viz.: 
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H™(A,/3) 



' dU?(\,!3) ^ _ 1 dU™(\,(3) _ 
d(3 Tx cos/3 d\ T(3 
1 dU?(\,/3) dU?(\,/3) 
cos/3 <9A A 9/3 p 



(4) 



The pair of vectors {T\,Tg} define the local tangential coordinate system in the plane or- 
thogonal to s , directed north and east, respectively. Spherical harmonics U™ are counted 
by degrees n — 0, 1, . . . and orders m = —n, —n + 1, . . . , n. Explicitly, 



where P™ are the associated Legendre polynomials. The first pair of vector harmonics are 
generated from the scalar harmonic £7°, with the electric component cos /3r^ and the magnetic 
component cos/3t\. The common index j used in Eq. [31 introduced for simplicity, counts 
all individual harmonics in the following manner: for each degree n all orders of electric 
harmonics are lined up, followed by all orders of magnetic harmonics. A particular vector 
harmonic, V 6 , which is the magnetic harmonic Hj -1 is depicted in Fig. [TJ It is equivalent to 
a left-handed rotation around the pole at A = 90°, /3 = 0. 



The most important advantage of the spherical harmonics is that they are orthogonal on 
the unit sphere in the space of continuous functions, or nearly orthogonal when discretized on 
a large set of uniformly distributed points. In the latter case, the deviation from orthogonality 
is negligibly small for a sufficiently large number of grid stars (> 10 4 ) of uniform density on 
the sky. The degree of uniformity and the number of grid stars should be higher if the higher 
order harmonics are to be used in the direct global solution. Normally, the lower orders of 
spherical harmonics are of interest, where the largest correlated errors emerge. Therefore, 
the series in Eq. [3] can be truncated in the new condition equations for the fitting coefficients 
Sj, Mj, and Qj. The solution for a subset of model terms is exact if the terms are orthogonal. 
In practice, the degree of orthogonality should be verified for a given distribution of stars and 
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2.2. A direct solution on a laptop 
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Fig. 1. — The magnetic vector spherical harmonic H x 
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weights (in a weighted LS). Rearranging the discretized spherical harmonics for the three 
types of astrometric unknowns as columns of a design matrix, the unknown coefficients and 
the right-hand side data as column vectors, the LS problem can be written in this compact 
matrix form: 

" S " 

[V s V M U ro ] M = y (6) 

. Q . 

The length of the design matrix is still very large in this setup, because it includes all the 
observations of grid stars. For JMAPS, it is about 144 million for 1 million grid stars. The 
width of the design matrix, on the other hand, is defined by how many spherical harmonic 
terms we want to solve for. Indeed, a sufficiently accurate solution can be obtained for 
any subset of spherical harmonics, as long as the columns of the design matrix are nearly 
orthogonal. This is verified by computing the correlation coefficients of the covariance matrix. 
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They should all be small, e.g., less than 0.01 in absolute value. If this is the case, including 
more terms in the design matrix will not change the solution for Sj, Mj, and Qj significantly. 
The number of unknowns can be made comfortably small for a given computer. We found, 
for example, that a global solution can be obtained within 1 hour for 400 unknowns on a 
regular laptop computer. 

Even with a limited number of unknowns, the design matrix is too large to be handled in 
fast memory without swapping with disk. However, there is no need to keep the entire matrix 
in memory if the observations are sorted by time. If the design matrix D = [V s V M U ro ] is 
divided into a number of blocks Dj in the vertical dimension (not necessarily of the same 
length), the normal matrix is the sum of the normal sub- matrices, D T D = Y2j DjDj- The 
accumulated normal matrix can be easily inverted due to its relatively small size, resulting in 
the covariance matrix, Cov = (D T D) . The off-diagonal elements of Cov should be small 
due to the near-orthogonality of the model terms, unless some additional global parameters 
are included. The diagonal elements are the variances of unit weight of the coefficients 
of spherical harmonics Sj, Mj, and Qj. If the observations are weighted by the expected 
standard deviation of measurement error, the variances are the squares of the standard 
errors carried by the corresponding spherical harmonics. The total mission-average variance 
is approximately the sum of the variances of the complete set of harmonics for each of the 
5 astrometric parameters. Since we obtain the variances for a limited set of harmonics, 
the total mission- average error can not be inferred from this computation. However, the 
uncertainty of specific harmonic components is accurately computed. 



3. Results and Discussion 

JMAPS is a pointing astrometric telescope with a single viewing direction. Without the 
ability of Hipparcos to simultaneously observe stars that are far apart on the sky, the re- 
quired rigidity of the reference system and the accuracy of astrometric parameters is achieved 
through measuring a number of carefully selected ICRF and radio-mute quasars and other 
extragalactic objects. These objects provide absolute constraints on positions (using the 
ICRF coordinates of superior accuracy), parallaxes and proper motions, which are negligibly 
small because of the extreme remoteness. The entire sky is observed with a 4-fold overlap. 
Astrometric observations are normally made around the great circle perpendicular to the 
sun direction. Some 72 observations per object are expected to be collected in three years. 
For the simulations described in this paper, we used a catalog of 44 ICRF quasars and 80 
compact extragalactic sources (QSO), which are not in the ICRF. All these reference quasars 
are brighter than magnitude 15. Only a subset of all observable stars, usually between 1 and 
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2 million strong, is used in the direct global solution. 

The attitude unknowns are represented by three parameters for each frame, viz., the 
translations along the axes of the detector and rotation around the boresight vector. These 
unknowns are eliminated frame by frame, reducing the number of conditions by three. The 
instrument calibration unknowns are represented in these simulations by sets of up to 28 
Zernike polynomials of field coordinates, separately for either coordinate in the detector 
plane. The first Zernike polynomials, which are constant, are excluded to avoid deficiency 
of rank, because they are indistinguishable from the attitude translations. In our simplified 
simulations, the calibration parameters are assumed to be constant within calibration blocks 
of equal length. Usually, blocks of 92 or 96 consecutive frames are used, corresponding to 
roughly 50 min of uninterrupted observations. As soon as a complete calibration block is 
collected, the QR factorization is applied, and the remaining astrometric equations are pre- 
multiplied with Q T , as well as the right-hand side. The number of condition equations is 
further reduced by 2 Nz for each calibration block, with Nz being the number of calibration 
parameters. This algorithm allows us to include a set of global parameters, which do not vary 
with time, such as the PPN 7-parameter. The accuracy or precision of global parameters can 
be reliably estimated, because they are mostly correlated with the low-order components. 
The number of vector spherical unknowns is 2(Ny + 1) 2 — 2, where Ny is the limiting degree, 
and the number of scalar spherical harmonics (for parallax) is {Ny + l) 2 . 



3.1. Accidental errors 



The statistical properties of accidental correlated errors are defined by the global co- 
variance matrix of the coefficients of spherical harmonics. The diagonal of the covariance 
matrix at Ny = 7 includes 126 vector spherical harmonic coefficients for positions and proper 
motions each and 64 scalar spherical harmonic coefficients for parallax. The square roots of 
the portions of the diagonal corresponding to each astrometric parameter are the standard 
deviations of error of unit weight, represented by a particular harmonic. For example, the 
standard deviation of the parallax zero-point error is the the standard deviation of the first 
spherical harmonic coefficient multiplied by the weighted average single measurement preci- 
sion of stars and quasars. Figs. [2] a and b show the standard deviations of harmonic errors 
of proper motions and parallax, respectively, obtained from a typical simulation of JMAPS 
mission. Generally, we find that the correlated errors in all three parameters fairly rapidly 
decline with the degree of spherical harmonic. To use an analogy from spectroscopy, in that 
sense, the spectrum of accidental errors is "red". There are some obvious "spectral lines", 
however, which are caused by the observing pattern and the distribution of reference quasars 



on the sky. 
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Fig. 2. — Standard deviations of the coefficients of spherical harmonics representing acciden- 
tal error in a) proper motions, and b) parallax. The errors are normalized to unit weight. 



We find that the distribution of accidental error becomes flatter with a significantly 
larger number of grid quasars. The relative height of the "spectral lines" depends on the 
distribution of grid quasars on the sky and the composition of the calibration model. Large 
holes in the distribution of quasars cause considerable degradation of the overall performance. 
Using the near-orthogonality of the discretized spherical function, the variance of accidental 
error of, e.g., proper motion at a given point (A,/3) can be estimated as 

var[MA,/3)] ~ J>ar [Mj] || V,(A, /3) || 2 = ^ C^\\ V^A, /3) f (7) 

3 3 

where is the corresponding part of the covariance matrix. Since this decomposition is 
limited to a finite set of spherical functions, only a lower bound of the total error can be 
obtained. Still, the distribution of the error carried by the lower-order harmonics is very in- 
formative. For example, one can estimate the degree of inhomogeneity of the correlated error 
on the celestial sphere, which can be significant for JMAPS. Fig. [3] depicts the distribution 
of the standard deviation of the total accidental error of parallax, which is contained in the 
first 64 scalar spherical harmonics. The plot is rotated into the Galactic coordinate system 
to emphasize the impact of the zone of avoidance around the Galactic plane where quasars 
brighter than 15 mag can not be found. The quasars, which were used to constrain the global 
solution for parallax, are shown as black dots. The build-up of error in the areas devoid of 
grid quasars is quite obvious. As a way of improving the overall mission performance, fainter 
quasars should be found near the Galactic plane and included in the grid. 
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Fig. 3. — Sky distribution of the standard deviation of parallax unit weight error for a 
simulated JMAPS mission. The reference quasars used to constrain the parallax solution 
are marked with black dots. 

It should be noted that the distribution in Fig. |3] corresponds to the expectancy of the 
spatially correlated error rather than to the outcome of a single mission. In other words, it 
shows what would emerge if the same mission is simulated many times with random noise 
in the input data, and the sample variance of the resulting errors is computed at each point. 
A single random realization of correlated error can be computed by 

6fi(\,0)= (c|r) T V(A,/3) (8) 

i 

with C/i being the unique positive definite matrix square root of C M , r a random vector drawn 
from Af(Q, 1/iVy) of N v elements, and V(A,/3) the column vector of the values Vj(\,f3). 
These transformations are performed separately for each coordinate direction for the vector- 
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valued parameters (position and proper motion). 

3.2. Systematic errors 

Systematic errors of global solutions are much harder to predict and analyze, because 
there are multiple sources of such errors, which are rarely known beforehand. Slowly varying 
perturbations of observational data, caused by external circumstances, are of special interest, 
as they can bring about smooth, large-scale errors. The orientation of the astrometric 
satellite with respect to the sun direction is one of the conceivable sources of systematic 
error. The angle between the sun direction, which is confined to the ecliptic plane, and 
the viewing direction changes in a predetermined way, because the entire celestial sphere 
should be observed as uniformly as possible. The thermal flow inevitably changes inside the 
telescope, resulting in slowly varying instrument parameters, e.g., the effective focal length 
or the basic angle for Gaia. If these variations are correlated with the celestial coordinates, 
there is no averaging out of the perturbation, and the error can propagate into the final 
catalog. In many cases, such specific physical influences can not be accurately modeled or 
predicted. A more general modeling approach can be exploited, where a certain perturbation 
is represented as a set of basis functions. For example, a systematic variation of the basic 
angle can be represented as a Fourier series of the sun angle, and each of the Fourier terms can 
be simulated separately. The previous studies for Hipparcos and SIM indicate that many of 
such elementary perturbations are benign, in that they cause a relatively small error. There 
are, however, some particularly dangerous perturbations, which may propagate into the final 
catalog with considerable magnification. Such harmful systematic effects should be identified 
and mitigated, if possible. This requires numerous mission simulations with different initial 
data, which may not be feasible for the extremely computer-intensive solutions for millions 
of individual grid objects. The proposed technique is fast enough to be used for massive 
simulations of slowly varying systematic perturbations, when the emerging astrometric error 
is confined to the lower degrees of spherical harmonics. 

Figure H] shows the results of a specific simulation for JMAPS, where a perturbation in 
the field-dependent calibration parameter Z2 was injected in the observational data, but not 
fitted out in the global solution. The term Z2 (second Zernike polynomial) corresponds to the 
differential scale of the instrument. The magnitude of the perturbation was normalized to 1 
mas at the edge of the field of view. The simulated observations (X and Y measurements) 
were free of random noise, to see more clearly the emerging pattern of the correlated error. 
The absolute error in the coefficients of 64 lower-order spherical harmonics is depicted in 
Fig. Hk. The spectrum of the error is dominated by the harmonic number 7 (which is U®), 
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Fig. 4. — Systematic error of JMAPS parallax resulting from a constant perturbation of the 
differential scale by 1 mas at the edge of the field of view: a) in the coefficients of the first 
64 spherical harmonics; b) on the sky. All values are in mas. 

followed by harmonic 19 (U^ 2 ) and so on. The total absolute error of parallax at a given 
point is the sum of all spherical harmonic errors. The total error in the first 64 harmonics is 
depicted in Fig. Sb. It shows that the simulated perturbation is one of the harmful errors for 
JMAPS, because it compounds to a perturbation of up to 2.6 mas in some parts of the sky, 
which is larger than the initial magnitude. Clearly, the distribution of constraining quasars, 
shown with black dots, plays a major role in the propagation of this systematic error, which 
compounds to larger values in the areas where the quasars are few. If the calibration term 
Z 2 is included in the set of fitting parameters in the global solution, the emerging error is 
zero in the absence of random noise. 



4. Conclusions 

We developed a method to investigate the properties of very large astrometric solutions, 
which involve unknown parameters for millions of celestial sources, as well as millions of 
nuisance unknowns. The method is based on a stepwise elimination of the attitude and 
calibration unknowns and the replacement of individual astrometric corrections with their 
expansions in orthogonal spherical functions of celestial coordinates. This approach works 
well for the JMAPS and SIM missions and could potentially be useful for Hipparcos and 
Gaia. However, demonstrating the applicability of the method to Hipparcos-like missions 
would require considerable adjustments, mostly related to the dynamic character of attitude 
parameters, which is beyond the scope of this paper. In particular, fitting a set of dynamic 



- 16 - 



parameters for each extended interval of uninterrupted rotation may render the proposed 
technique of QR-elimination of the attitude unknowns impractical. An additional complica- 
tion arises for Gaia, where each of the multiple CCDs in the focal plane requires a separate set 
of calibration parameters. The pointing, or step-stare mode of operation of JMAPS makes 
it best-suited for the proposed global solution technique with block-wise elimination of at- 
titude and calibration parameters, so that complete analysis for realistic sky coverages and 
observing schedules can be performed for billions of condition equations. The propagation 
of zonal and correlated errors of both accidental and systematic origin can be successfully 
computed using this method. When the number of expansion terms is appropriately small, 
full mission solutions can be obtained using regular computers within a few hours with more 
than a hundred million unknowns. Some applications of the spherical harmonic solution to 
the JMAPS mission are described. 
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